#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "../include/var_struc.h"
#include "../include/finite_volume.h"


int fluid_var_update(struct flu_var * FV, struct cell_var * cv)
{
	const int dim = (int)config[0];		
	const int num_cell = (int)config[3];
	struct i_f_var ifv;

	for(int k = 0; k < num_cell; ++k)
		{
			cons_qty_copy_cv2ifv(&ifv, cv, k);			
			if(cons2prim(&ifv) == 0)
				{
					fprintf(stderr, "\n Wrong in copying cons_var to prim_var!\n");
					fprintf(stderr, "ifv.prim: %d,\tRHO=%.10lf,\tP=%.10lf,\tPHI=%.10lf,\tZ_a=%.10lf\n",k,ifv.RHO,ifv.P,ifv.PHI,ifv.Z_a);
					fprintf(stderr, "ifv.cons: %d,\tU_rho=%.10lf,\tU_e=%.10lf,\tU_phi=%.10lf,\tU_u=%.10lf,\tU_v=%.10lf\n",k,ifv.U_rho,ifv.U_e,ifv.U_phi,ifv.U_u,ifv.U_v);
					return 0;
				}
			cv->U_e_a[k] = ifv.U_e_a; 
			prim_var_copy_ifv2FV(&ifv, FV, k);

			cons_qty_copy_ifv2cv(&ifv, cv, k);			
		}
	return 1;
}


static int order2_i_f_var_init(const struct cell_var * cv, struct i_f_var * ifv, const int k)
{
	const int dim = (int)config[0];	
	const double n_x = ifv->n_x, n_y = ifv->n_y, n_z = ifv->n_z;
	const double delta_x = ifv->delta_x, delta_y = ifv->delta_y, delta_z = ifv->delta_z;
	const double eps = config[4];

	ifv->d_rho  = cv->gradx_rho[k]*n_x;
	ifv->d_e    = cv->gradx_e[k]  *n_x;
	ifv->d_u    = cv->gradx_u[k]  *n_x;
	if ((int)config[2] == 2)
		{					
			ifv->d_phi  = cv->gradx_phi[k]*n_x;
			ifv->d_z_a  = cv->gradx_z_a[k]*n_x;
		}
	if (dim > 1)
		{
			ifv->d_rho += cv->grady_rho[k]*n_y;
			ifv->d_e   += cv->grady_e[k]  *n_y;
			ifv->d_u   += cv->grady_u[k]  *n_y;
			ifv->d_v    = cv->gradx_v[k]  *n_x + cv->grady_v[k]*n_y;
			if ((int)config[2] == 2)
				{
					ifv->d_phi += cv->grady_phi[k]*n_y;
					ifv->d_z_a += cv->grady_z_a[k]*n_y;
				}
		}
	if (dim > 2)
		{
			ifv->d_rho += cv->gradz_rho[k]*n_z;
			ifv->d_e   += cv->gradz_e[k]  *n_z;
			ifv->d_u   += cv->gradz_u[k]  *n_z;
			ifv->d_v   += cv->gradz_v[k]  *n_z;
			ifv->d_w    = cv->gradx_w[k]  *n_x + cv->grady_w[k]*n_y + cv->gradz_w[k]*n_z;
			if ((int)config[2] == 2)
				{
					ifv->d_phi += cv->gradz_phi[k]*n_z;
					ifv->d_z_a += cv->gradz_z_a[k]*n_z;					
				}
		}

	ifv->t_rho  = cv->gradx_rho[k]*(-n_y);
	ifv->t_e    = cv->gradx_e[k]  *(-n_y);
	ifv->t_u    = cv->gradx_u[k]  *(-n_y);
	if ((int)config[2] == 2)
		{					
			ifv->t_phi  = cv->gradx_phi[k]*(-n_y);
			ifv->t_z_a  = cv->gradx_z_a[k]*(-n_y);
		}
	if (dim > 1)
		{
			ifv->t_rho += cv->grady_rho[k]*n_x;
			ifv->t_e   += cv->grady_e[k]  *n_x;
			ifv->t_u   += cv->grady_u[k]  *n_x;
			ifv->t_v    = cv->gradx_v[k]  *(-n_y) + cv->grady_v[k]*n_x;
			if ((int)config[2] == 2)
				{
					ifv->t_phi += cv->grady_phi[k]*n_x;
					ifv->t_z_a += cv->grady_z_a[k]*n_x;
				}
		}

	if (cons2prim(ifv) == 0)
		{
			fprintf(stderr, "Error happens on primitive variable!\n");
			return 0;
		}
	if ((int)config[31] == 0)
		{					
			ifv->d_p = ifv->d_e;
			ifv->t_p = ifv->t_e;
			
			ifv->RHO += cv->gradx_rho[k]*delta_x;
			ifv->P   += cv->gradx_e[k]  *delta_x;
			ifv->U   += cv->gradx_u[k]  *delta_x;
			if ((int)config[2] == 2)
				{									
					ifv->PHI += cv->gradx_phi[k]*delta_x;
					ifv->Z_a += cv->gradx_z_a[k]*delta_x;
				}
			if (dim > 1)
				{
					ifv->RHO += cv->grady_rho[k]*delta_y;
					ifv->P   += cv->grady_e[k]  *delta_y;
					ifv->U   += cv->grady_u[k]  *delta_y;
					ifv->V   += cv->gradx_v[k]  *delta_x + cv->grady_v[k]*delta_y;
					if ((int)config[2] == 2)
						{
							ifv->PHI += cv->grady_phi[k]*delta_y;							
							ifv->Z_a += cv->grady_z_a[k]*delta_y;
						}				
				}
			if (dim > 2)
				{
					ifv->RHO += cv->gradz_rho[k]*delta_z;
					ifv->P   += cv->gradz_e[k]  *delta_z;
					ifv->U   += cv->gradz_u[k]  *delta_z;
					ifv->V   += cv->gradz_v[k]  *delta_z;
					ifv->W   += cv->gradx_w[k]  *delta_x + cv->grady_w[k]*delta_y + cv->gradz_w[k]*delta_z;
					if ((int)config[2] == 2)
						{
							ifv->PHI += cv->gradz_phi[k]*delta_z;				
							ifv->Z_a += cv->gradz_z_a[k]*delta_z;							
						}				
				}
			ifv->gamma = 1.0+1.0/(ifv->Z_a/(config[6]-1.0)+(1.0-ifv->Z_a)/(config[106]-1.0));
			if (ifv->P < eps || ifv->RHO < eps || ifv->PHI < -eps || ifv->PHI > (1.0+eps) || ifv->Z_a < -eps || ifv->Z_a > (1.0+eps))
			{
				printf("\n Out of variable range after reconstruction!\n");
				printf("ifv->P=%lf, ifv->RHO=%lf, ifv->PHI=%.10lf, ifv->Z_a=%.10lf\n", ifv->P, ifv->RHO, ifv->PHI, ifv->Z_a);
			}
		}
	else if ((int)config[31] == 1)
		{
			ifv->d_u = (ifv->d_u - ifv->U*ifv->d_rho)/ifv->RHO;
			ifv->d_p = (ifv->d_e - 0.5*ifv->d_rho*ifv->U*ifv->U - ifv->RHO*ifv->U*ifv->d_u) * (ifv->gamma-1.0);	
			if (dim > 1)
				{
					ifv->d_v  = (ifv->d_v - ifv->V*ifv->d_rho)/ifv->RHO;
					ifv->d_p += (- 0.5*ifv->d_rho*ifv->V*ifv->V - ifv->RHO*ifv->V*ifv->d_v) * (ifv->gamma-1.0);
				}
			if (dim > 2)
				{			
					ifv->d_w  = (ifv->d_w - ifv->W*ifv->d_rho)/ifv->RHO;
					ifv->d_p += (- 0.5*ifv->d_rho*ifv->W*ifv->W - ifv->RHO*ifv->W*ifv->d_w) * (ifv->gamma-1.0);
				}
			if ((int)config[2] == 2)
				ifv->d_phi = (ifv->d_phi - ifv->PHI*ifv->d_rho)/ifv->RHO;
			if (!isinf(config[60]))
				ifv->d_gamma = (ifv->d_gamma - ifv->gamma*ifv->d_rho)/ifv->RHO;

			ifv->U_rho += cv->gradx_rho[k]*delta_x;
			ifv->U_e   += cv->gradx_e[k]  *delta_x;
			ifv->U_u   += cv->gradx_u[k]  *delta_x;
			if ((int)config[2] == 2)									
				ifv->U_phi += cv->gradx_phi[k]*delta_x;
			if (dim > 1)
				{
					ifv->U_rho += cv->grady_rho[k]*delta_y;
					ifv->U_e   += cv->grady_e[k]  *delta_y;
					ifv->U_u   += cv->grady_u[k]  *delta_y;
					ifv->U_v   += cv->gradx_v[k]  *delta_x + cv->grady_v[k]*delta_y;
					if ((int)config[2] == 2)				
						ifv->U_phi += cv->grady_phi[k]*delta_y;				
				}
			if (dim > 2)
				{
					ifv->U_rho += cv->gradz_rho[k]*delta_z;
					ifv->U_e   += cv->gradz_e[k]  *delta_z;
					ifv->U_u   += cv->gradz_u[k]  *delta_z;
					ifv->U_v   += cv->gradz_v[k]  *delta_z;
					ifv->U_w   += cv->gradx_w[k]  *delta_x + cv->grady_w[k]*delta_y + cv->gradz_w[k]*delta_z;
					if ((int)config[2] == 2)
						ifv->U_phi += cv->gradz_phi[k]*delta_z;				
				}
			if(cons2prim(ifv) == 0)
				{
					fprintf(stderr, "Error happens on primitive variable!\n");
					return 0;
				}
		}
	return 1;
}


static int order2_i_f_var0(struct i_f_var * ifv)
{
	const int dim = (int)config[0];
		
	ifv->d_rho = 0.0;
	ifv->d_e   = 0.0;
	ifv->d_u   = 0.0;
	if (dim > 1)						
		ifv->d_v = 0.0;
	if (dim > 2)
		ifv->d_w = 0.0;	

	ifv->t_rho = 0.0;
	ifv->t_e   = 0.0;
	ifv->t_u   = 0.0;
	ifv->t_v = 0.0;
			
	if ((int)config[2] == 2)
		{
			ifv->d_phi = 0.0;
			ifv->d_z_a = 0.0;
			ifv->t_phi = 0.0;
			ifv->t_z_a = 0.0;
		}

	if(cons2prim(ifv) == 0)
		{
			fprintf(stderr, "Error happens on primitive variable!\n");
			return 0;
		}
	
	return 1;
}

	  
int interface_var_init(const struct cell_var * cv, const struct mesh_var * mv,
					   struct i_f_var * ifv, struct i_f_var * ifv_R,
					   const int k, const int j, const int i)
{
	const int dim = (int)config[0];
	const int order = (int)config[9];
	int **cc = cv->cell_cell;
	int **cp = mv->cell_pt;

	int p_p, p_n;
	if(j == cp[k][0]-1) 
		{
			p_p=cp[k][1];
			p_n=cp[k][j+1];
		}				  
	else
		{
			p_p=cp[k][j+2];
			p_n=cp[k][j+1];
		}

	if (dim == 1)
		ifv->n_x = 1.0;
	if (dim > 1)
		{
			ifv->n_x = cv->n_x[k][j];
			ifv->n_y = cv->n_y[k][j];
		}
	if (dim > 2)
		ifv->n_z = cv->n_z[k][j];
	if (dim == 2)
		ifv->length = sqrt((mv->X[p_p] - mv->X[p_n])*(mv->X[p_p] - mv->X[p_n]) + (mv->Y[p_p] - mv->Y[p_n])*(mv->Y[p_p] - mv->Y[p_n]));

	cons_qty_copy_cv2ifv(ifv, cv, k);

	if (order == 2)
		{
			if (dim == 1)
				ifv->delta_x = mv->X[p_n] - cv->X_c[k];
			else if (dim == 2)
				{					
					ifv->delta_x = 0.5*(mv->X[p_p] + mv->X[p_n]) - cv->X_c[k];
					ifv->delta_y = 0.5*(mv->Y[p_p] + mv->Y[p_n]) - cv->Y_c[k];
				}			
			if(order2_i_f_var_init(cv, ifv, k) == 0)			
				{
					fprintf(stderr, "Error happens on primitive variable!\n");
					return 0;
				}
		}
		
	ifv_R->n_x = ifv->n_x;
	if (dim > 1)
		{
			ifv_R->n_y = ifv->n_y;
			ifv_R->length = ifv->length;
		}
	if (dim > 2)
		ifv_R->n_z = ifv->n_z;
	
	int cR; //cell_right	
	if (cc[k][j] >= 0)
		{
			cR = cc[k][j];
			cons_qty_copy_cv2ifv(ifv_R, cv, cR);			

			if (order == 2)
				{
					if (dim == 1)
						ifv_R->delta_x = mv->X[p_n] - cv->X_c[cR];
					else if (dim == 2)
						{
							ifv_R->delta_x = 0.5*(mv->X[p_p] + mv->X[p_n]) - cv->X_c[cR];
							ifv_R->delta_y = 0.5*(mv->Y[p_p] + mv->Y[p_n]) - cv->Y_c[cR];
						}
					if(order2_i_f_var_init(cv, ifv_R, cR) == 0)
						{
							fprintf(stderr, "Error happens on primitive variable!\n");
							return 0;
						}
				}
		}
	else if (cc[k][j] == -1)//initial boundary condition.
		{
			if (i > 0)
				return -1;
			cons_qty_copy_cv2ifv(ifv_R, cv, k);

			if (order == 2)
				if(order2_i_f_var0(ifv_R) == 0)
					{
						fprintf(stderr, "Error happens on primitive variable!\n");
						return 0;
					}
		}
	else if (cc[k][j] == -3)//prescribed boundary condition.
		{
			cons_qty_copy_cv2ifv(ifv_R, cv, k);			

			if (order == 2)
				if(order2_i_f_var0(ifv_R) == 0)
					{
						fprintf(stderr, "Error happens on primitive variable!\n");
						return 0;
					}
		}		
	else if (cc[k][j] != -2&&cc[k][j] != -4)
		{
			printf("No suitable boundary!cc = %d!\n",cc[k][j]);
			return 0;
		}

	if (order == 1)
		{
			if(cons2prim(ifv) == 0)
				{
					fprintf(stderr, "Error happens on primitive variable!\n");
					return 0;
				}
			if (cc[k][j] != -2&&cc[k][j] != -4)
				if(cons2prim(ifv_R) == 0)
					{
						fprintf(stderr, "Error happens on primitive variable!\n");
						return 0;
					}
		}

	double u_R, v_R;
	if (cc[k][j] == -2)//reflecting boundary condition.
		{
			*ifv_R = *ifv;
			u_R =  ifv_R->U*ifv_R->n_x + ifv_R->V*ifv_R->n_y;
			u_R = -u_R;
			v_R = -ifv_R->U*ifv_R->n_y + ifv_R->V*ifv_R->n_x;
			ifv_R->U = u_R*ifv_R->n_x - v_R*ifv_R->n_y;
			ifv_R->V = u_R*ifv_R->n_y + v_R*ifv_R->n_x;
		}
	else if (cc[k][j] == -4)//symmetry boundary condition.
		*ifv_R = *ifv;

	return 1;
}


double tau_calc(const struct cell_var * cv, const struct mesh_var * mv)
{
	const double CFL = config[7];
	if (CFL < 0.0)
		return -CFL;
	const int dim = (int)config[0];		
	const int num_cell = (int)config[3];
	int ** cp = mv->cell_pt;
	
	double tau = config[1];
	struct i_f_var ifv, ifv_R;
	double cum, lambda_max;
	int ivi;
	
	double qn, qn_R;
	double c, c_R;	
	
	for(int k = 0; k < num_cell; ++k)
		{
			cum = 0.0;
			
			for(int j = 0; j < cp[k][0]; ++j)
				{
					ivi = interface_var_init(cv, mv, &ifv, &ifv_R, k, j, 0);
					if (ivi < 0)
						;
					else if(ivi == 0)
						return -1.0;
					else if (dim == 1)
						{
							c = sqrt(ifv.gamma * ifv.P / ifv.RHO);
							c_R = sqrt(ifv_R.gamma * ifv_R.P / ifv_R.RHO);
							lambda_max = fmax(c+fabs(ifv.U), c_R+fabs(ifv_R.U));
							cum += lambda_max;	
						}
					else if (dim == 2)
						{
							qn = ifv.U*ifv.n_x + ifv.V*ifv.n_y; 
							qn_R = ifv_R.U*ifv_R.n_x + ifv_R.V*ifv_R.n_y;
							c = sqrt(ifv.gamma * ifv.P / ifv.RHO);
							c_R = sqrt(ifv_R.gamma * ifv_R.P / ifv_R.RHO);
							lambda_max = fmax(c+fabs(qn), c_R+fabs(qn_R));
							cum += 0.5*lambda_max * ifv.length;
						}
				}
			tau = fmin(tau, cv->vol[k]/cum * CFL);
		}	//To decide tau.
	return tau;
}

#define CV_MOVE(var)  cv->var[i] = cv->var[j]
#define CV_MOVE2(var)				\
    do {					\
	cv->var[i][0] = cv->var[j][0];		\
	cv->var[i][1] = cv->var[j][1];		\
	cv->var[i][2] = cv->var[j][2];		\
	cv->var[i][3] = cv->var[j][3];		\
    } while (0)

void mesh_offset(const struct cell_var * cv, const struct mesh_var * mv)
{	
    int i;
    for(i = 0; i < mv->num_pt; i++)
	{
	    mv->X[i] += config[210];
	    mv->Y[i] += config[211];
	}
	const int num_cell = (int)config[3];
	const int num_cell_ghost = mv->num_ghost + num_cell;
    for(i = 0; i < num_cell_ghost; i++)
	{
	    cv->X_c[i] += config[210];
	    cv->Y_c[i] += config[211];
	}
}
	
void mesh_move(const struct cell_var * cv, const struct mesh_var * mv)
{	
	int i, j;
	for(i = 0; i < mv->num_pt; i++)
	    mv->X[i] += config[10];
	const int num_cell = (int)config[3];
	const int num_cell_ghost = mv->num_ghost + num_cell;  // Sod_mesh (mv.num_ghost = 0)
	for(i = 0; i < num_cell_ghost; i++)
		cv->X_c[i] += config[10];

	const int dim = (int)config[0];
	const int n_x = (int)config[13], n_y = (int)config[14];
	for(i = 0; i < num_cell; i++)
	    {
		if (i%n_x == n_x-1)
		    j = i;
		else
		    j = i+1;
		CV_MOVE(U_rho);			
		CV_MOVE(U_gamma);
		CV_MOVE(U_e);
		CV_MOVE(U_u);
		CV_MOVE2(RHO_p);
		CV_MOVE2(U_p);
		CV_MOVE2(P_p);
		if (dim > 1)
		    {
			CV_MOVE(U_v);
			CV_MOVE2(V_p);
		    }
		if (dim > 2)
		    CV_MOVE(U_w);
		if ((int)config[2] == 2)
		    {
			CV_MOVE(U_phi);
			CV_MOVE(U_e_a);
			CV_MOVE2(PHI_p);
			CV_MOVE2(Z_a_p);
		    }
	    }
}
